clear all;
small = 1.0e-10;

% -- File Directories   
figdir = 'fig/';
outdir = 'out/';
matdir = 'mat/';
mcmcdir = '../matlab/mcmc/';

% Model Name
mod_name = '_baseline';


% Load Posterior Elements
fstr = [mcmcdir 'F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'm_draws' mod_name];load(fstr);
fstr = [mcmcdir 'g_matrix' mod_name];load(fstr);
fstr = [mcmcdir 'sigma2_F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'mtau_F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'stddev_del50_m_draws' mod_name];load(fstr);
fstr = [mcmcdir 'F_rho_draws' mod_name];load(fstr);
fstr = [mcmcdir 'mu_U_draws' mod_name];load(fstr);

sig_delta_a_draws = sqrt(sigma2_F_draws(2,:));
sig_m_draws = sqrt(sigma2_F_draws(1,:));
hl_draws = log(0.5)./log(F_rho_draws);
trnd_g = g_matrix(:,2);
dt_g = trnd_g(2)-trnd_g(1); % Scaling used in UM's program
m_mean_draws = mtau_F_draws(2,:)*dt_g;

n_t = 118;
calvec = (1900:1:2017)';
n_f = size(F_draws,1);
n_rep = size(F_draws,2);

% Load fitted trends (for plots)
fstr = [matdir 'yfit_mat'];load(fstr);
yfit_mat_countries = yfit_mat(:,1:end-1);
n_c = size(yfit_mat_countries,2);

load([matdir 'country_code']);
load([matdir 'oecd_list']);
i_oecd = zeros(n_c,1);
nn = size(oecd_list,1);
for j = 1:nn;
    ss = char(oecd_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_oecd(jj) = 1;
    end;
end;

% Fitted Value of f(t) .. low-frequency projection
Fproj_draws = g_matrix*F_draws;
Fproj_pct = NaN(n_t,3);
pct = [1/6 0.50 5/6];
for t = 1:n_t;
    tmp = Fproj_draws(t,:)';
    tmp_p = pctile(tmp,pct);
    Fproj_pct(t,:) = tmp_p';
end;

% m(t) draws
mproj_draws = g_matrix*m_draws; % This is m(t) without its mean
mproj_draws = dif(mproj_draws,1);
mproj_draws = mproj_draws + repmat(m_mean_draws,n_t,1);
mproj_pct = NaN(n_t,3);
pct = [1/6 0.50 5/6];
for t = 1:n_t;
    tmp = mproj_draws(t,:)';
    tmp_p = pctile(tmp,pct);
    mproj_pct(t,:) = tmp_p';
end;

% Figure 3_a_b
figure_3_a_b;

% Figure 3_c;
sig_m_min = 0.001;
sig_m_max = 0.02;
sig_m_grid = linspace(sig_m_min,sig_m_max,25)';
sig_m_prior = NaN(25,1);
del = 1/(13*13);
tmp = del*(1:12)';
sig_m_prior(1:12) = tmp;
sig_m_prior(13) = del*13;
sig_m_prior(14:25) = flipud(tmp);
% Get posterior probabilities
sig_m_post = NaN(size(sig_m_grid));
for i = 1:size(sig_m_grid,1);
    sig_m_post(i) = mean( abs(sig_m_draws' - sig_m_grid(i)) < .0001);
end;

figure_3_c;

% Results for Table 2_a
pct = [0.05 1/6 0.50 5/6 0.95];
sig_delta_a_pct = 100*pctile(sig_delta_a_draws',pct)';
mu_m_pct = 100*pctile(m_mean_draws',pct)';
sig_m_pct = 100*pctile(sig_m_draws',pct)';
hl_pct = pctile(hl_draws',pct)';
mu_c_pct = pctile(mu_U_draws',pct)';

outfile_name = [outdir 'Table_2_a.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Parameter, ');
prtmat_comma(pct,fileID,'%5.2f','\n');
fprintf(fileID,'Sig_delta_a,');
prtmat_comma(sig_delta_a_pct,fileID,'%5.2f','\n');
fprintf(fileID,'mu_m,');
prtmat_comma(mu_m_pct,fileID,'%5.2f','\n');
fprintf(fileID,'Sig_m,');
prtmat_comma(sig_m_pct,fileID,'%5.2f','\n');
fprintf(fileID,'hl,');
prtmat_comma(hl_pct,fileID,'%5.0f','\n');
fprintf(fileID,'\n mu_c,');
prtmat_comma(mu_c_pct,fileID,'%5.2f','\n');

% Write f quantiles CSV File

outfile_name = [outdir 'F_Insample.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'In-sample posterior quantiles for F \n');
fprintf(fileID,'Year,,1/6, 0.5, 5/6 \n');
for t = 1:n_t;
    fprintf(fileID,'%4i ,, ',calvec(t));
    prtmat_comma(Fproj_pct(t,:),fileID,'%5.3f','\n');
end;






